Ejemplo 6

Inicialización

Datos globales para modelos

#--- Ejemplo 6 ---
#-Reading data-
desastres<-read.table("desastres.txt",header=TRUE)
n<-nrow(desastres)
plot(desastres,type="l")

plot(desastres[2:n,2]-desastres[1:(n-1),2],type="l")

plot(log(desastres[2:n,2])-log(desastres[1:(n-1),2]),type="l")

#-Defining data-
data<-list("n"=n,"y"=desastres$No.Desastres,"x"=desastres$Anho)
data<-list("n"=n,"y"=c(desastres$No.Desastres[1:(n-6)],rep(NA,6)),"x"=desastres$Anho)

Funciones de graficación

source("util.R")
graficarY1 <- function(modelo){
  if("BUGSoutput" %in% names(modelo)){
    modelo <- modelo$BUGSoutput
  }
  out.sum<-modelo$summary
  #Predictions
  out.yf<-out.sum[grep("yf1",rownames(out.sum)),]
  ymin<-min(desastres[,2],out.yf[,c(1,3,7)])
  ymax<-max(desastres[,2],out.yf[,c(1,3,7)])
  
  par(mfrow=c(1,1))
  plot(desastres,type="l",col="grey80",ylim=c(ymin,ymax))
  lines(desastres[,1],out.yf[,1],lwd=2,col=2)
  lines(desastres[,1],out.yf[,3],lty=2,col=2)
  lines(desastres[,1],out.yf[,7],lty=2,col=2)
  lines(desastres[,1],out.yf[,5],lwd=2,col=4)
  
  #Medias
  out.mu<-out.sum[grep("mu",rownames(out.sum)),]
  par(mfrow=c(1,1))
  plot(desastres,type="l",col="grey80")
  lines(desastres[,1],out.mu[,1],lwd=2,col=2)
  lines(desastres[,1],out.mu[,3],lty=2,col=2)
  lines(desastres[,1],out.mu[,7],lty=2,col=2)
}
graficarTasa <- function(modelo){
  if("BUGSoutput" %in% names(modelo)){
    modelo <- modelo$BUGSoutput
  }
  out.sum<-modelo$summary
  out.tasa<-out.sum[grep("lambda",rownames(out.sum)),]
  out.tasa<-out.sum[grep("p",rownames(out.sum)),]
  or<-order(mortality$x)
  ymin<-min(mortality$y/mortality$n,out.tasa[,c(1,3,7)])
  ymax<-max(mortality$y/mortality$n,out.tasa[,c(1,3,7)])
  plot(mortality$x,mortality$y/mortality$n,ylim=c(ymin,ymax))
  lines(mortality$x[or],out.tasa[or,1],lwd=2,col=4)
  lines(mortality$x[or],out.tasa[or,3],lty=2,col=4)
  lines(mortality$x[or],out.tasa[or,7],lty=2,col=4)
}

Ejercicio 6a

cat Ej6a.txt
## model
## {
##   #Likelihood
##   for (i in 1:n) {
##     #Neg Binomial
##     #    y[i] ~ dnegbin(p[i],r)
##     #    logit(p[i])<-beta[1]+beta[2]*x[i]
##     #    mu[i]<-r*(1-p[i])/p[i]
##     #Poisson
##      y[i] ~ dpois(mu[i])
##      log(mu[i])<-beta[1]+beta[2]*x[i]
##   }
##   #Priors 
##   for (j in 1:2) { beta[j] ~ dnorm(0,0.001) }
##   #Neg Binomial
##   #aux ~ dpois(5)
##   #r <- aux + 1
##   
##   #Prediction 1
##   #Neg Binomial
##   #for (i in 1:n) { yf1[i] ~ dnegbin(p[i],r) }
##   #Poisson
##   for (i in 1:n) { yf1[i] ~ dpois(mu[i]) }
## 
## }
#-Defining inits-
inits<-function(){list(beta=rep(0,2),yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux=1,aux2=1,yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux2=1,yf1=rep(1,n),tau.y=1)}
# inits<-function(){list(beta=rep(0,n),tau.b=1,yf1=rep(1,n))}
# inits<-function(){list(mu=rep(1,n),tau.b=1,yf1=rep(1,n))}

#-Selecting parameters to monitor-
parameters<-c("beta","yf1","mu")
# parameters<-c("beta","yf1","mu","tau")
# parameters<-c("beta","yf1","mu","tau","tau.y")
# parameters<-c("beta","yf1","mu","r")
# parameters<-c("beta","yf1","mu","tau","r")
# parameters<-c("tau.b","yf1","mu")
ej6a.bugs<-bugs(data,inits,parameters,model.file="Ej6a.txt",
               n.iter=20000,n.chains=1,n.burnin=4000)
ej6a.jags<-jags(data,inits,parameters,model.file="Ej6a.txt",
              n.iter=20000,n.chains=4,n.burnin=4000,n.thin=5)
## module glm loaded
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 106
##    Unobserved stochastic nodes: 120
##    Total graph size: 679
## 
## Initializing model
### es mejor poner un adelgazamiento cuando hay alta correlación entre \Beta_1 y \Beta_2
print("Bugs")
## [1] "Bugs"
exploracionMCMC(ej6a.bugs)

print("Jags")
## [1] "Jags"
exploracionMCMC(ej6a.jags)

Betas

betasMCMC(ej6a.bugs)
##                mean     2.5%    97.5% prob
## beta[1] 34.97886250 25.40000 44.99000    0
## beta[2] -0.01816264 -0.02347 -0.01311    0
betasMCMC(ej6a.jags)
##                mean        2.5%      97.5% prob
## beta[1] 35.04765163 25.49703245 44.9492355    0
## beta[2] -0.01819738 -0.02343192 -0.0131291    0

DIC

dicMCMC(ej6a.bugs)
## [1] 338.1
dicMCMC(ej6a.jags)
## [1] 339.0566
graficarY1(ej6a.bugs)

graficarY1(ej6a.jags)

Ejercicio 6a - Predicción

cat Ej6a.txt
## model
## {
##   #Likelihood
##   for (i in 1:n) {
##     #Neg Binomial
##     #    y[i] ~ dnegbin(p[i],r)
##     #    logit(p[i])<-beta[1]+beta[2]*x[i]
##     #    mu[i]<-r*(1-p[i])/p[i]
##     #Poisson
##      y[i] ~ dpois(mu[i])
##      log(mu[i])<-beta[1]+beta[2]*x[i]
##   }
##   #Priors 
##   for (j in 1:2) { beta[j] ~ dnorm(0,0.001) }
##   #Neg Binomial
##   #aux ~ dpois(5)
##   #r <- aux + 1
##   
##   #Prediction 1
##   #Neg Binomial
##   #for (i in 1:n) { yf1[i] ~ dnegbin(p[i],r) }
##   #Poisson
##   for (i in 1:n) { yf1[i] ~ dpois(mu[i]) }
## 
## }
data<-list("n"=n,"y"=c(desastres$No.Desastres[1:(n-6)],rep(NA,6)),"x"=desastres$Anho)
#-Defining inits-
inits<-function(){list(beta=rep(0,2),yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux=1,aux2=1,yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux2=1,yf1=rep(1,n),tau.y=1)}
# inits<-function(){list(beta=rep(0,n),tau.b=1,yf1=rep(1,n))}
# inits<-function(){list(mu=rep(1,n),tau.b=1,yf1=rep(1,n))}

#-Selecting parameters to monitor-
parameters<-c("beta","yf1","mu")
# parameters<-c("beta","yf1","mu","tau")
# parameters<-c("beta","yf1","mu","tau","tau.y")
# parameters<-c("beta","yf1","mu","r")
# parameters<-c("beta","yf1","mu","tau","r")
# parameters<-c("tau.b","yf1","mu")
ej6a.bugs<-bugs(data,inits,parameters,model.file="Ej6a.txt",
               n.iter=20000,n.chains=1,n.burnin=4000)
ej6a.jags<-jags(data,inits,parameters,model.file="Ej6a.txt",
              n.iter=20000,n.chains=4,n.burnin=4000,n.thin=5)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 106
##    Unobserved stochastic nodes: 120
##    Total graph size: 679
## 
## Initializing model
### es mejor poner un adelgazamiento cuando hay alta correlación entre \Beta_1 y \Beta_2
print("Bugs")
## [1] "Bugs"
exploracionMCMC(ej6a.bugs)

print("Jags")
## [1] "Jags"
exploracionMCMC(ej6a.jags)

Betas

betasMCMC(ej6a.bugs)
##                mean     2.5%    97.5% prob
## beta[1] 34.97886250 25.40000 44.99000    0
## beta[2] -0.01816264 -0.02347 -0.01311    0
betasMCMC(ej6a.jags)
##                mean        2.5%       97.5%       prob
## beta[1] 35.05034708 25.64597116 45.11129741 7.8125e-05
## beta[2] -0.01819942 -0.02351621 -0.01321869 7.8125e-05

DIC

dicMCMC(ej6a.bugs)
## [1] 338.1
dicMCMC(ej6a.jags)
## [1] 339.3104

GrƔficas

graficarY1(ej6a.bugs)

graficarY1(ej6a.jags)

Restore data

#reestablecer para no tener predicción
data<-list("n"=n,"y"=desastres$No.Desastres,"x"=desastres$Anho)

Ejercicio 6aa

cat Ej6aa.txt
## model
## {
##   #Likelihood
##   for (i in 1:n) {
##     #Neg Binomial
##      y[i] ~ dnegbin(p[i],r)
##      logit(p[i])<-beta[1]+beta[2]*x[i]
##      mu[i]<-r*(1-p[i])/p[i]
##     #Poisson
##     #    y[i] ~ dpois(mu[i])
##     #    log(mu[i])<-beta[1]+beta[2]*x[i]
##   }
##   #Priors 
##   for (j in 1:2) { beta[j] ~ dnorm(0,0.001) }
##   #Neg Binomial
##   aux ~ dpois(10)
##   r <- aux + 1
##   
##   #Prediction 1
##   #Neg Binomial
##   for (i in 1:n) { yf1[i] ~ dnegbin(p[i],r) }
##   #Poisson
##   #for (i in 1:n) { yf1[i] ~ dpois(mu[i]) }
## 
## }
#-Defining inits-
inits<-function(){list(beta=rep(0,2),yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux=1,aux2=1,yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux2=1,yf1=rep(1,n),tau.y=1)}
# inits<-function(){list(beta=rep(0,n),tau.b=1,yf1=rep(1,n))}
# inits<-function(){list(mu=rep(1,n),tau.b=1,yf1=rep(1,n))}

#-Selecting parameters to monitor-
# parameters<-c("beta","yf1","mu")
# parameters<-c("beta","yf1","mu","tau")
# parameters<-c("beta","yf1","mu","tau","tau.y")
parameters<-c("beta","yf1","mu","r")
# parameters<-c("beta","yf1","mu","tau","r")
# parameters<-c("tau.b","yf1","mu")
ej6a.bugs<-bugs(data,inits,parameters,model.file="Ej6aa.txt",
               n.iter=10000,n.chains=2,n.burnin=1000)
ej6a.jags<-jags(data,inits,parameters,model.file="Ej6aa.txt",
              n.iter=10000,n.chains=2,n.burnin=1000,n.thin=1)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 112
##    Unobserved stochastic nodes: 115
##    Total graph size: 1244
## 
## Initializing model
### es mejor poner un adelgazamiento cuando hay alta correlación entre \Beta_1 y \Beta_2

GrƔficas

print("Bugs")
## [1] "Bugs"
exploracionMCMC(ej6a.bugs)

print("Jags")
## [1] "Jags"
exploracionMCMC(ej6a.jags)

exploracionMCMC_r(ej6a.bugs)

exploracionMCMC_r(ej6a.jags)

Betas

betasMCMC(ej6a.bugs)
##                 mean      2.5%    97.5% prob
## beta[1] -33.72750056 -43.80000 -23.9200    0
## beta[2]   0.01878268   0.01356   0.0241    0
betasMCMC(ej6a.jags)
##                 mean         2.5%       97.5% prob
## beta[1] -18.04765706 -25.63400634 -5.63592620    0
## beta[2]   0.01047011   0.00391693  0.01451908    0

DIC

dicMCMC(ej6a.bugs)
## [1] 347.5122
dicMCMC(ej6a.jags)
## [1] 399.8575
graficarY1(ej6a.bugs)

graficarY1(ej6a.jags)

Ejercicio 6b

cat Ej6b.txt
## model
## {
## #Likelihood
## for (i in 1:n) {
## #Neg Binomial
## #    y[i] ~ dnegbin(p[i],r)
## #    logit(p[i])<-beta[1]+beta[2]*step(x[i]-tau)
## #    mu[i]<-r*(1-p[i])/p[i]
## #Poisson
##  y[i] ~ dpois(mu[i])
##  log(mu[i])<-beta[1]+beta[2]*step(x[i]-tau)
##  }
## #Priors 
## for (j in 1:2) { beta[j] ~ dnorm(0,0.001) }
## #Neg Binomial
## #aux ~ dpois(5)
## #r <- aux + 1
## aux2 ~ dcat(a[])
## tau <- aux2 + 1850
## for (j in 1:112) { a[j]<- 1/112}
## #Prediction 1
## #Neg Binomial
## #for (i in 1:n) { yf1[i] ~ dnegbin(p[i],r) }
## #Poisson
## for (i in 1:n) { yf1[i] ~ dpois(mu[i]) }
## 
## }
#-Defining inits-
inits<-function(){list(beta=rep(0,2),yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux=1,aux2=1,yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux2=1,yf1=rep(1,n),tau.y=1)}
# inits<-function(){list(beta=rep(0,n),tau.b=1,yf1=rep(1,n))}
# inits<-function(){list(mu=rep(1,n),tau.b=1,yf1=rep(1,n))}

#-Selecting parameters to monitor-
parameters<-c("beta","yf1","mu")
# parameters<-c("beta","yf1","mu","tau")
# parameters<-c("beta","yf1","mu","tau","tau.y")
# parameters<-c("beta","yf1","mu","r")
# parameters<-c("beta","yf1","mu","tau","r")
# parameters<-c("tau.b","yf1","mu")
ej6b.bugs<-bugs(data,inits,parameters,model.file="Ej6b.txt",
               n.iter=5000,n.chains=1,n.burnin=500)
ej6b.jags<-jags(data,inits,parameters,model.file="Ej6b.txt",
              n.iter=5000,n.chains=1,n.burnin=500,n.thin=5)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 112
##    Unobserved stochastic nodes: 115
##    Total graph size: 1133
## 
## Initializing model
### es mejor poner un adelgazamiento cuando hay alta correlación entre \Beta_1 y \Beta_2
print("Bugs")
## [1] "Bugs"
exploracionMCMC(ej6b.bugs)

print("Jags")
## [1] "Jags"
exploracionMCMC(ej6b.jags)

Betas

betasMCMC(ej6b.bugs)
##              mean      2.5%     97.5% prob
## beta[1]  1.142430  0.951585  1.320000    0
## beta[2] -1.261346 -1.570525 -0.965045    0
betasMCMC(ej6b.jags)
##              mean       2.5%      97.5%        prob
## beta[1]  1.078548  0.9452293  1.3186034 0.005555556
## beta[2] -1.170457 -1.5511344 -0.9268255 0.005555556

DIC

dicMCMC(ej6b.bugs)
## [1] 341.1391
dicMCMC(ej6b.jags)
## [1] 366.9928
graficarY1(ej6b.bugs)

graficarY1(ej6b.jags)

Ejercicio 6bb

cat Ej6bb.txt
## model
## {
## #Likelihood
## for (i in 1:n) {
## #Neg Binomial
##  y[i] ~ dnegbin(p[i],r)
##  logit(p[i])<-beta[1]+beta[2]*step(x[i]-tau)
##  mu[i]<-r*(1-p[i])/p[i]
## #Poisson
## #    y[i] ~ dpois(mu[i])
## #    log(mu[i])<-beta[1]+beta[2]*step(x[i]-tau)
##  }
## #Priors 
## for (j in 1:2) { beta[j] ~ dnorm(0,0.001) }
## #Neg Binomial
## aux ~ dpois(5)
## r <- aux + 1
## aux2 ~ dcat(a[])
## tau <- aux2 + 1850
## for (j in 1:112) { a[j]<- 1/112}
## #Prediction 1
## #Neg Binomial
## for (i in 1:n) { yf1[i] ~ dnegbin(p[i],r) }
## #Poisson
## #for (i in 1:n) { yf1[i] ~ dpois(mu[i]) }
## 
## }
#-Defining inits-
inits<-function(){list(beta=rep(0,2),yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux=1,aux2=1,yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux2=1,yf1=rep(1,n),tau.y=1)}
# inits<-function(){list(beta=rep(0,n),tau.b=1,yf1=rep(1,n))}
# inits<-function(){list(mu=rep(1,n),tau.b=1,yf1=rep(1,n))}

#-Selecting parameters to monitor-
parameters<-c("beta","yf1","mu")
# parameters<-c("beta","yf1","mu","tau")
# parameters<-c("beta","yf1","mu","tau","tau.y")
# parameters<-c("beta","yf1","mu","r")
# parameters<-c("beta","yf1","mu","tau","r")
# parameters<-c("tau.b","yf1","mu")
ej6b.bugs<-bugs(data,inits,parameters,model.file="Ej6bb.txt",
               n.iter=5000,n.chains=2,n.burnin=1000)
ej6b.jags<-jags(data,inits,parameters,model.file="Ej6bb.txt",
              n.iter=5000,n.chains=2,n.burnin=1000,n.thin=5)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 112
##    Unobserved stochastic nodes: 116
##    Total graph size: 1810
## 
## Initializing model
### es mejor poner un adelgazamiento cuando hay alta correlación entre \Beta_1 y \Beta_2
print("Bugs")
## [1] "Bugs"
exploracionMCMC(ej6b.bugs)

print("Jags")
## [1] "Jags"
exploracionMCMC(ej6b.jags)

Betas

betasMCMC(ej6b.bugs)
##              mean      2.5%    97.5%    prob
## beta[1] 0.8680623 0.2378725 1.415025 0.00425
## beta[2] 1.2579773 0.9203925 1.590000 0.00000
betasMCMC(ej6b.jags)
##              mean      2.5%    97.5%    prob
## beta[1] 0.8469051 0.2165423 1.427110 0.00375
## beta[2] 1.2590502 0.9339628 1.603908 0.00000

DIC

dicMCMC(ej6b.bugs)
## [1] 345.9282
dicMCMC(ej6b.jags)
## [1] 346.0538
graficarY1(ej6b.bugs)

graficarY1(ej6b.jags)

Ejercicio 6c

cat Ej6c.txt
## model
## {
## #Likelihood
## #Space eq.
## for (i in 1:n) {
##  y[i] ~ dpois(mu[i])
##  log(mu[i])<-beta[i]
##  }
## #State eq.
## for (i in 2:n) {
##  beta[i] ~ dnorm(beta[i-1],tau.b)
##  }
## #Priors 
## beta[1] ~ dnorm(0,0.001)
## tau.b ~ dgamma(0.001,0.001)
## 
## #Prediction 1
## for (i in 1:n) { yf1[i] ~ dpois(mu[i]) }
## 
## }
#-Defining inits-
# inits<-function(){list(beta=rep(0,2),yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux=1,aux2=1,yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux2=1,yf1=rep(1,n),tau.y=1)}
inits<-function(){list(beta=rep(0,n),tau.b=1,yf1=rep(1,n))}
# inits<-function(){list(mu=rep(1,n),tau.b=1,yf1=rep(1,n))}

#-Selecting parameters to monitor-
# parameters<-c("beta","yf1","mu")
# parameters<-c("beta","yf1","mu","tau")
# parameters<-c("beta","yf1","mu","tau","tau.y")
# parameters<-c("beta","yf1","mu","r")
# parameters<-c("beta","yf1","mu","tau","r")
parameters<-c("beta", "tau.b","yf1","mu")
ej6c.bugs<-bugs(data,inits,parameters,model.file="Ej6c.txt",
               n.iter=5000,n.chains=1,n.burnin=500)
ej6c.jags<-jags(data,inits,parameters,model.file="Ej6c.txt",
              n.iter=5000,n.chains=1,n.burnin=500,n.thin=5)
## Warning in jags.model(model.file, data = data, inits = init.values,
## n.chains = n.chains, : Unused variable "x" in data
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 112
##    Unobserved stochastic nodes: 225
##    Total graph size: 454
## 
## Initializing model
### es mejor poner un adelgazamiento cuando hay alta correlación entre \Beta_1 y \Beta_2
print("Bugs")
## [1] "Bugs"
exploracionMCMC(ej6c.bugs)

print("Jags")
## [1] "Jags"
exploracionMCMC(ej6c.jags)

Betas

betasMCMC2 <- function(modelo){
  if("BUGSoutput" %in% names(modelo)){
    modelo <- modelo$BUGSoutput
  }
  out <- modelo$sims.list
  out.sum<-modelo$summary
  out.sum.t<-out.sum[grep("beta",rownames(out.sum)),c(1,3,7)]
  out.sum.t<-cbind(out.sum.t,apply(out$beta,2,prob))
  dimnames(out.sum.t)[[2]][4]<-"prob"
  # print(modelo$sims.list$beta)
  print(out.sum.t)
}
betasMCMC2(ej6c.bugs)
##                   mean        2.5%       97.5%         prob
## beta[1]    1.186404476  0.64511750  1.71700000 0.0004444444
## beta[2]    1.176569556  0.67658493  1.66252492 0.0000000000
## beta[3]    1.124161800  0.64859499  1.57500000 0.0000000000
## beta[4]    1.049636000  0.56629499  1.48052492 0.0000000000
## beta[5]    1.019035578  0.54708000  1.44000000 0.0002222222
## beta[6]    1.055533889  0.58808997  1.48352492 0.0000000000
## beta[7]    1.066987567  0.62878493  1.48852492 0.0000000000
## beta[8]    1.077790578  0.63284248  1.49000000 0.0000000000
## beta[9]    1.064095044  0.60964750  1.47505000 0.0002222222
## beta[10]   1.120860822  0.68908997  1.54100000 0.0000000000
## beta[11]   1.114888267  0.68958493  1.53700000 0.0000000000
## beta[12]   1.111782956  0.69667482  1.49452492 0.0000000000
## beta[13]   1.112582609  0.68458997  1.50500000 0.0000000000
## beta[14]   1.087464256  0.66088250  1.48052500 0.0004444444
## beta[15]   1.134025000  0.71293746  1.52152492 0.0000000000
## beta[16]   1.206312933  0.79754249  1.61300000 0.0000000000
## beta[17]   1.213023956  0.80384249  1.62500000 0.0000000000
## beta[18]   1.230625867  0.83744428  1.62704969 0.0000000000
## beta[19]   1.261067356  0.86422233  1.66252492 0.0000000000
## beta[20]   1.256854356  0.87554750  1.67852493 0.0000000000
## beta[21]   1.242330578  0.85782738  1.65852492 0.0000000000
## beta[22]   1.197264644  0.79614249  1.59000000 0.0000000000
## beta[23]   1.162509911  0.75122230  1.54752492 0.0000000000
## beta[24]   1.178227378  0.76949499  1.57900000 0.0000000000
## beta[25]   1.179947733  0.77742231  1.56900000 0.0000000000
## beta[26]   1.164413933  0.75654750  1.55100000 0.0000000000
## beta[27]   1.203783556  0.79613746  1.61852492 0.0000000000
## beta[28]   1.206855889  0.79754937  1.62652492 0.0000000000
## beta[29]   1.172528311  0.76814750  1.59652492 0.0000000000
## beta[30]   1.145873956  0.73849499  1.58000000 0.0000000000
## beta[31]   1.097180844  0.67338493  1.49852492 0.0000000000
## beta[32]   1.078060978  0.66419499  1.48757425 0.0000000000
## beta[33]   1.004883378  0.57760000  1.41452491 0.0000000000
## beta[34]   0.950060960  0.51099499  1.36552491 0.0000000000
## beta[35]   0.911446511  0.46928500  1.34352500 0.0002222222
## beta[36]   0.865600564  0.41949500  1.30705000 0.0002222222
## beta[37]   0.781966251  0.33398000  1.22000000 0.0013333333
## beta[38]   0.700140861  0.23969000  1.11905000 0.0015555556
## beta[39]   0.641424687  0.18442250  1.08452500 0.0037777778
## beta[40]   0.557892641  0.06231775  0.99670500 0.0157777778
## beta[41]   0.467437037 -0.05962300  0.92395250 0.0377777778
## beta[42]   0.393509395 -0.17106750  0.86500750 0.0711111111
## beta[43]   0.332014646 -0.22832000  0.80730500 0.1057777778
## beta[44]   0.279232468 -0.27471500  0.77456000 0.1437777778
## beta[45]   0.233980967 -0.32697750  0.73136750 0.1811111111
## beta[46]   0.197431397 -0.35985250  0.69696250 0.2191111111
## beta[47]   0.118350283 -0.44276250  0.61862500 0.3093333333
## beta[48]   0.063620082 -0.52201250  0.57589750 0.3786666667
## beta[49]   0.037567481 -0.54387750  0.55365250 0.4197777778
## beta[50]   0.007925561 -0.57476250  0.52636750 0.4682222222
## beta[51]  -0.005325450 -0.61120500  0.52350000 0.4820000000
## beta[52]  -0.017300148 -0.63160000  0.53476750 0.4902222222
## beta[53]  -0.021949274 -0.63070500  0.53286250 0.4733333333
## beta[54]  -0.002374363 -0.57310500  0.55223000 0.4964444444
## beta[55]   0.039074907 -0.55871000  0.61164500 0.4484444444
## beta[56]   0.031045862 -0.54865000  0.59510000 0.4577777778
## beta[57]   0.032573767 -0.53162000  0.59450000 0.4413333333
## beta[58]   0.060416980 -0.48085250  0.61166000 0.4064444444
## beta[59]   0.041604198 -0.51738250  0.57896750 0.4340000000
## beta[60]   0.003886435 -0.55991000  0.54470000 0.4828888889
## beta[61]  -0.055997085 -0.61332500  0.46666250 0.4322222222
## beta[62]  -0.093295581 -0.64683500  0.42205750 0.3753333333
## beta[63]  -0.134151626 -0.69679250  0.37190500 0.3017777778
## beta[64]  -0.176386264 -0.75173500  0.34095250 0.2568888889
## beta[65]  -0.222535966 -0.79571000  0.29851000 0.2191111111
## beta[66]  -0.250927162 -0.85855250  0.27646750 0.1902222222
## beta[67]  -0.283007907 -0.89337250  0.26708750 0.1617777778
## beta[68]  -0.295794442 -0.89886750  0.25710500 0.1375555556
## beta[69]  -0.313690998 -0.91796250  0.19672000 0.1228888889
## beta[70]  -0.310900799 -0.91280500  0.20661000 0.1246666667
## beta[71]  -0.284486115 -0.87571000  0.21991000 0.1457777778
## beta[72]  -0.243625697 -0.82836250  0.24407250 0.1822222222
## beta[73]  -0.235973155 -0.85269750  0.25327750 0.1908888889
## beta[74]  -0.233601357 -0.85687500  0.26327250 0.1900000000
## beta[75]  -0.216706930 -0.87686750  0.26910500 0.2146666667
## beta[76]  -0.180229579 -0.79441500  0.29446250 0.2617777778
## beta[77]  -0.125199487 -0.71162750  0.37233250 0.3275555556
## beta[78]  -0.073217687 -0.63300500  0.45077750 0.3915555556
## beta[79]  -0.025089270 -0.56238750  0.51972000 0.4540000000
## beta[80]   0.055059698 -0.46825250  0.61167250 0.4246666667
## beta[81]   0.113543070 -0.39249750  0.66156250 0.3391111111
## beta[82]   0.153783501 -0.34525250  0.70487750 0.2891111111
## beta[83]   0.147676497 -0.34041000  0.69402500 0.3040000000
## beta[84]   0.150433047 -0.34150250  0.69286250 0.3068888889
## beta[85]   0.155751608 -0.34341500  0.70704500 0.2977777778
## beta[86]   0.146421960 -0.36305750  0.70730500 0.3237777778
## beta[87]   0.141963327 -0.35497250  0.69420500 0.3168888889
## beta[88]   0.145182792 -0.33694000  0.68863500 0.3162222222
## beta[89]   0.154931453 -0.33313500  0.70361000 0.3066666667
## beta[90]   0.168015509 -0.35680500  0.73011000 0.2915555556
## beta[91]   0.162448044 -0.35905750  0.74795000 0.3004444444
## beta[92]   0.088263030 -0.44461500  0.64138750 0.3833333333
## beta[93]  -0.012043187 -0.56198250  0.52694500 0.4815555556
## beta[94]  -0.089814005 -0.66776250  0.45510500 0.3864444444
## beta[95]  -0.147530418 -0.74976250  0.39850000 0.3208888889
## beta[96]  -0.181250019 -0.82110000  0.37100750 0.2728888889
## beta[97]  -0.213382314 -0.85020000  0.33967250 0.2315555556
## beta[98]  -0.324215623 -1.02105000  0.23230000 0.1326666667
## beta[99]  -0.419107433 -1.15657500  0.16766250 0.0835555556
## beta[100] -0.499654283 -1.30157500  0.10521000 0.0533333333
## beta[101] -0.558980694 -1.40857500  0.04982850 0.0404444444
## beta[102] -0.630100846 -1.47552500 -0.00975905 0.0226666667
## beta[103] -0.685981070 -1.55352500 -0.03760975 0.0184444444
## beta[104] -0.733273178 -1.64052500 -0.05942825 0.0137777778
## beta[105] -0.765173364 -1.67500000 -0.11404250 0.0084444444
## beta[106] -0.785009173 -1.76757500 -0.13663750 0.0062222222
## beta[107] -0.795944328 -1.80200000 -0.13360000 0.0051111111
## beta[108] -0.822018586 -1.87605000 -0.15024750 0.0060000000
## beta[109] -0.845735502 -1.92062500 -0.13918750 0.0051111111
## beta[110] -0.856390510 -1.98262500 -0.10969750 0.0077777778
## beta[111] -0.878551151 -2.08205000 -0.10609000 0.0091111111
## beta[112] -0.888979448 -2.11862500 -0.08120125 0.0142222222
betasMCMC2(ej6c.jags)
##                    mean        2.5%        97.5%        prob
## beta[1]    1.163482e+00  0.61653332  1.688036414 0.000000000
## beta[2]    1.155103e+00  0.69062415  1.661144926 0.000000000
## beta[3]    1.104094e+00  0.62047084  1.542039619 0.000000000
## beta[4]    1.036272e+00  0.55949126  1.486278335 0.000000000
## beta[5]    1.012751e+00  0.49804356  1.459594189 0.001111111
## beta[6]    1.059754e+00  0.62721248  1.472797977 0.000000000
## beta[7]    1.084537e+00  0.69323320  1.508018336 0.000000000
## beta[8]    1.098809e+00  0.69791268  1.486678246 0.000000000
## beta[9]    1.086438e+00  0.67627994  1.512372381 0.000000000
## beta[10]   1.144797e+00  0.75348161  1.535781203 0.000000000
## beta[11]   1.129157e+00  0.73143817  1.525132094 0.000000000
## beta[12]   1.120300e+00  0.70079123  1.504650873 0.000000000
## beta[13]   1.114881e+00  0.70143074  1.497270777 0.000000000
## beta[14]   1.073313e+00  0.64792592  1.438048293 0.000000000
## beta[15]   1.117922e+00  0.69650869  1.481535467 0.000000000
## beta[16]   1.189446e+00  0.79861095  1.590036396 0.000000000
## beta[17]   1.185427e+00  0.73011400  1.598368200 0.000000000
## beta[18]   1.203445e+00  0.81508297  1.585723633 0.000000000
## beta[19]   1.240088e+00  0.81907015  1.697939469 0.000000000
## beta[20]   1.235901e+00  0.82777764  1.679552441 0.000000000
## beta[21]   1.217214e+00  0.79915733  1.616699970 0.000000000
## beta[22]   1.171920e+00  0.78091496  1.533623247 0.000000000
## beta[23]   1.136336e+00  0.71349773  1.536792642 0.000000000
## beta[24]   1.161959e+00  0.75188490  1.575042077 0.000000000
## beta[25]   1.167846e+00  0.76108877  1.574735262 0.000000000
## beta[26]   1.150074e+00  0.75437353  1.545563593 0.000000000
## beta[27]   1.192170e+00  0.81183412  1.596889877 0.000000000
## beta[28]   1.191645e+00  0.83814811  1.588207551 0.000000000
## beta[29]   1.164381e+00  0.76911834  1.592200963 0.000000000
## beta[30]   1.129464e+00  0.77283331  1.528607882 0.000000000
## beta[31]   1.080048e+00  0.70840429  1.464420587 0.000000000
## beta[32]   1.051264e+00  0.66552333  1.457800381 0.000000000
## beta[33]   9.813956e-01  0.59429529  1.396188831 0.000000000
## beta[34]   9.266961e-01  0.51515735  1.329929522 0.000000000
## beta[35]   8.829496e-01  0.44839284  1.303133477 0.000000000
## beta[36]   8.359225e-01  0.41095324  1.294551129 0.000000000
## beta[37]   7.492307e-01  0.33146463  1.157445192 0.000000000
## beta[38]   6.740897e-01  0.23423120  1.110356486 0.002222222
## beta[39]   6.145917e-01  0.17257305  1.021404194 0.003333333
## beta[40]   5.361867e-01  0.05319936  0.943473359 0.013333333
## beta[41]   4.498741e-01 -0.03710060  0.875984208 0.042222222
## beta[42]   3.819056e-01 -0.15427053  0.802354605 0.073333333
## beta[43]   3.259579e-01 -0.18514443  0.795971909 0.088888889
## beta[44]   2.815060e-01 -0.21435915  0.724862365 0.110000000
## beta[45]   2.330995e-01 -0.26295307  0.657488049 0.152222222
## beta[46]   1.921749e-01 -0.30005829  0.630028311 0.217777778
## beta[47]   1.136508e-01 -0.39280931  0.566258514 0.307777778
## beta[48]   5.962312e-02 -0.48034487  0.533535380 0.394444444
## beta[49]   3.093576e-02 -0.47426470  0.501936724 0.423333333
## beta[50]  -6.078716e-04 -0.57446119  0.485795194 0.478888889
## beta[51]  -4.092496e-03 -0.60872509  0.483078815 0.491111111
## beta[52]  -2.026932e-03 -0.56205485  0.508664677 0.477777778
## beta[53]   4.010552e-05 -0.59312281  0.521204908 0.475555556
## beta[54]   2.270509e-02 -0.56487598  0.528294140 0.465555556
## beta[55]   7.644045e-02 -0.45451180  0.572506719 0.365555556
## beta[56]   7.692947e-02 -0.43879035  0.567063112 0.368888889
## beta[57]   8.484146e-02 -0.44327232  0.565981430 0.348888889
## beta[58]   1.180576e-01 -0.41768240  0.578534875 0.294444444
## beta[59]   9.917300e-02 -0.44412396  0.589258202 0.310000000
## beta[60]   4.869441e-02 -0.50816172  0.563041565 0.386666667
## beta[61]  -1.430205e-02 -0.62021154  0.497161969 0.472222222
## beta[62]  -5.177988e-02 -0.65583266  0.438087809 0.461111111
## beta[63]  -1.038017e-01 -0.71483654  0.414743254 0.367777778
## beta[64]  -1.500083e-01 -0.74871144  0.351465317 0.301111111
## beta[65]  -2.065042e-01 -0.78662445  0.338876210 0.236666667
## beta[66]  -2.299535e-01 -0.86039571  0.298419469 0.222222222
## beta[67]  -2.790901e-01 -0.93307976  0.262859338 0.196666667
## beta[68]  -2.938798e-01 -0.96122647  0.254392398 0.184444444
## beta[69]  -3.174843e-01 -0.97663053  0.246526410 0.150000000
## beta[70]  -3.231845e-01 -1.01894244  0.248645459 0.152222222
## beta[71]  -3.074042e-01 -0.93658145  0.231949022 0.151111111
## beta[72]  -2.691884e-01 -0.90162595  0.267300871 0.190000000
## beta[73]  -2.629401e-01 -0.90989779  0.260624238 0.185555556
## beta[74]  -2.656030e-01 -0.92635537  0.232333012 0.190000000
## beta[75]  -2.535037e-01 -0.90792325  0.249824727 0.197777778
## beta[76]  -2.155179e-01 -0.82136378  0.291537114 0.244444444
## beta[77]  -1.542421e-01 -0.76890759  0.323221445 0.313333333
## beta[78]  -9.736057e-02 -0.68082126  0.401125747 0.377777778
## beta[79]  -4.160857e-02 -0.57936692  0.467998791 0.452222222
## beta[80]   3.918525e-02 -0.47499596  0.585720105 0.427777778
## beta[81]   1.018864e-01 -0.45137656  0.639037576 0.355555556
## beta[82]   1.453024e-01 -0.39149125  0.714086048 0.297777778
## beta[83]   1.382575e-01 -0.42858932  0.694016493 0.308888889
## beta[84]   1.306509e-01 -0.36596773  0.686906374 0.333333333
## beta[85]   1.433223e-01 -0.35410625  0.684472438 0.307777778
## beta[86]   1.383640e-01 -0.35133408  0.683252307 0.314444444
## beta[87]   1.375278e-01 -0.36178998  0.698959365 0.310000000
## beta[88]   1.437146e-01 -0.38701425  0.729229515 0.298888889
## beta[89]   1.511331e-01 -0.37243699  0.713111338 0.291111111
## beta[90]   1.754862e-01 -0.33596744  0.758854117 0.270000000
## beta[91]   1.732897e-01 -0.34962951  0.730623203 0.288888889
## beta[92]   1.142001e-01 -0.43139141  0.676031217 0.355555556
## beta[93]   1.366668e-02 -0.49220838  0.567634432 0.476666667
## beta[94]  -6.193843e-02 -0.58572179  0.464664670 0.417777778
## beta[95]  -1.060191e-01 -0.62126391  0.387466557 0.354444444
## beta[96]  -1.321277e-01 -0.67082218  0.360501039 0.340000000
## beta[97]  -1.555074e-01 -0.71298970  0.385579371 0.297777778
## beta[98]  -2.642365e-01 -0.83566182  0.247318327 0.201111111
## beta[99]  -3.475842e-01 -0.89264096  0.177701067 0.118888889
## beta[100] -4.208126e-01 -1.04027357  0.144216541 0.076666667
## beta[101] -4.899116e-01 -1.14609123  0.067911356 0.050000000
## beta[102] -5.653700e-01 -1.25525846  0.041526380 0.031111111
## beta[103] -6.329116e-01 -1.36135924 -0.005417412 0.024444444
## beta[104] -6.894248e-01 -1.44018748 -0.051513452 0.016666667
## beta[105] -7.259733e-01 -1.49806783 -0.059065080 0.011111111
## beta[106] -7.467424e-01 -1.44497269 -0.044509144 0.013333333
## beta[107] -7.537286e-01 -1.47905222 -0.037769138 0.016666667
## beta[108] -7.792557e-01 -1.55581617 -0.018457969 0.020000000
## beta[109] -8.001735e-01 -1.55813086 -0.054858819 0.014444444
## beta[110] -8.113561e-01 -1.62778758 -0.029250561 0.016666667
## beta[111] -8.337856e-01 -1.64641481 -0.049182394 0.016666667
## beta[112] -8.442536e-01 -1.81183085 -0.006933733 0.024444444
hist(ej6c.bugs$sims.list$beta)

DIC

dicMCMC(ej6c.bugs)
## [1] 338
dicMCMC(ej6c.jags)
## [1] 346.4019
graficarY1(ej6c.bugs)

graficarY1(ej6c.jags)

Ejercicio 6d

cat Ej6d.txt
## model
## {
## #Likelihood
## #Space eq.
## for (i in 1:n) {
##  y[i] ~ dpois(mu[i])
##  }
## #State eq.
## for (i in 2:n) {
##  mu[i] ~ dgamma(tau.b,b[i])
##  b[i] <- tau.b/mu[i-1]
##  }
## #Priors 
## mu[1] ~ dgamma(0.001,0.001)
## #tau.b ~ dgamma(0.001,0.001)
## tau.b ~ dgamma(200,1)
## 
## #Prediction 1
## for (i in 1:n) { yf1[i] ~ dpois(mu[i]) }
## 
## }
#-Defining inits-
# inits<-function(){list(beta=rep(0,2),yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux=1,aux2=1,yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,2),aux2=1,yf1=rep(1,n),tau.y=1)}
# inits<-function(){list(beta=rep(0,n),tau.b=1,yf1=rep(1,n))}
inits<-function(){list(mu=rep(1,n),tau.b=1,yf1=rep(1,n))}

#-Selecting parameters to monitor-
# parameters<-c("beta","yf1","mu")
# parameters<-c("beta","yf1","mu","tau")
# parameters<-c("beta","yf1","mu","tau","tau.y")
# parameters<-c("beta","yf1","mu","r")
# parameters<-c("beta","yf1","mu","tau","r")
parameters<-c("tau.b","yf1","mu")
ej6d.bugs<-bugs(data,inits,parameters,model.file="Ej6d.txt",
               n.iter=5000,n.chains=1,n.burnin=500)
ej6d.jags<-jags(data,inits,parameters,model.file="Ej6d.txt",
              n.iter=5000,n.chains=1,n.burnin=500,n.thin=5)
## Warning in jags.model(model.file, data = data, inits = init.values,
## n.chains = n.chains, : Unused variable "x" in data
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 112
##    Unobserved stochastic nodes: 225
##    Total graph size: 453
## 
## Initializing model
### es mejor poner un adelgazamiento cuando hay alta correlación entre \Beta_1 y \Beta_2
# print("Bugs")
# exploracionMCMC(ej6d.bugs)
# print("Jags")
# exploracionMCMC(ej6d.jags)

Betas

# betasMCMC(ej6d.bugs)
# betasMCMC(ej6d.jags)

DIC

dicMCMC(ej6d.bugs)
## [1] 343.8262
dicMCMC(ej6d.jags)
## [1] 342.861
# graficarY1(ej6d.bugs)
# graficarY1(ej6d.jags)

DICs

dicMCMC(ej6a.bugs)
## [1] 347.5122
dicMCMC(ej6a.jags)
## [1] 399.8575
dicMCMC(ej6b.bugs)
## [1] 345.9282
dicMCMC(ej6b.jags)
## [1] 346.0538
dicMCMC(ej6c.bugs)
## [1] 338
dicMCMC(ej6c.jags)
## [1] 346.4019
dicMCMC(ej6d.bugs)
## [1] 343.8262
dicMCMC(ej6d.jags)
## [1] 342.861